Supercell technique for total-energy calculations of finite charged and polar systems 
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We study the behaviour of total-energy supercell calculations for dipolar molecules and charged clus- 
ters. Using a cutoff Coulomb interaction within the framework of a plane-wave basis set formalism, 
with all other aspects of the method (pseudopotentials, basis set, exchange-correlation functional) 
unchanged, we are able to assess directly the interaction effects present in the supercell technique. 
We find that the supercell method gives structures and energies in almost total agreement with the 
results of calculations for finite systems, even for molecules with large dipole moments. We also show 
that the performance of finite-grid calculations can be improved by allowing a degree of aliasing in 
the Hartree energy, and by using a reciprocal space definition of the cutoff Coulomb interaction. 
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I. INTRODUCTION 

Density- functional theory (DFT), typically performed 
within a local approximation for exchange and corre- 
lation such as the local density approximation (LDA), 
has proved to be an extremely powerful and successful 
method for the calculation of ground-state propierties of 
electronic systems. The Kohn-Sham formalismlil allows 
the ground-state electronic density and total energy to be 
found through minimisation with respect to a set of one- 
electron orbitals. Expansion of these orbitals in terms 
of a plane-wave bp-sis set, in conjunction with the use of 
pseudopotentialsD to represent nuclei and core electrons, 
provides a systematic route to convergence for such calcu- 
lations and also allows the use of Fast Fourier Transforms 
(FFTs). Plane waves also make the calculation of forces 
for molecular dynamics straightforward. 

While the use of a plane- wave basis set is most natural 
for infinite periodic systems such as bulk solids, the com- 
putational advantages this offers apply equally for cal- 
culations of finite systems (molecules) or systems which 
are not infinite in all three dimensions (polymers or sur- 
faces). However, as the use of a discrete set of plane 
waves corresponds to employing periodic boundary con- 
ditions (PBCs), the calculation will still be essentially 
for an infinite system, with the supercell, containing the 
molecule, periodically repeated. These images will see 
each other through the long-ranged Coulomb interaction, 
thus introducing differences from the isolated system. 

In the case of a charged finite system, the supercell ap- 
proach gives a divergent Coulomb repulsion energy and 
so this interaction must be removed in some way, or alter- 
natively a fully finite calculation undertaken. A neutral 
system with a dipole (or higher order) moment will also 
have contributions to the total energy from the interac- 



tion between supercells: these will fall off with supercell 
size, but there is only power-law convergence in the en- 
ergy, and filling cells with extra vacuum may be imprac- 
tically expensive. Makov and PayneB developed separate 
corrections to the total energy functional for the cases of 
charged and dipolar systems, which improved the con- 
vergence of energy differences for test systems. 

In order to avoid the enforced periodicity of a supercell 
approach, many calculations for molecules and clusters 
are performed on a finite grid. In the case of a plane- 
wave basis this requires the use of computationally ex- 
pensive double-sized FFTs to perform the finite convo- 
lution for the Hartree potential. Alternatively, localized 
Gaussian basis sets are extensively employed in quan- 
tum chemistry, but these are harder to implement and 
do not allow for such straightforward systematic conver- 
gence procedures as plane waves. 

Finite-size effects are also relevant in other types of 
many-body calculations. Modified Coulomb interaction 
techniques have been used, for G W many-body pertur- 
bation theory calcuLatjonsQ and Quantum Monte Carlo 
(QMC) calculations,l3i3 and these authors proposed that 
a similar approach could be used for DFT studies of finite 
systems. However, the physical consequences of using a 
finite cell to approximate an infinite one differ, as the 
finite-size effects take the form of spurious correlations 
in QMC or screening from the contents of other cells in 
GW^work, and thus it is by no means clear that the rela- 
tive importance of removing the electrostatic interaction 
of the ground-state densities in DFT total-energy calcu- 
lations of finite systems will be as large. 

In this work we study charged and dipolar systems 
using a Hamiltonian with a cutoff Coulomb interaction, 
and compare with the results of a conventional supercell 
calculation. As we implement our modification for finite 
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systems within the framework of a pseudopotential plane- 
wave total-energy method (in this case CASTEPlI) , we are 
able to assess the relative performance, and specifically 
convergence, of these approaches directly, and investigate 
in which circumstances a fully isolated calculation is re- 
quired. Using the test cases of a Mg^ ion and a NaCl 
bond we demonstrate that our method removes all inter- 
action between supercells. We then study the cases of a 
Na2^ cluster and the polar molecules dimethylsulphox- 
ide and para-nitroanisole, which provide some idea of the 
error involved in unmoderated supercell calculations for 
realistic charged and polar finite systems. We find that 
there is almost total agreement in the ground-state geom- 
etry determined by the supercell method and for the iso- 
lated system, even for the smallest supercell sizes which 
can contain the density, though it has been suggested in 
the literatureB that supercell calculations under these cir- 
cumstances require substantial extra vacuum for conver- 
gence. All calculations were performed within the LDA 
(although our conclusions apply for any local approxi- 
mation to exchange and correlation), and sampled the 
Brillouin Zone at the F-point only. We use atomic units 
throughout. 



II. COMPUTATIONAL METHOD 

In this section we describe how to remove any interac- 
tion between supercells in a plane-wave basis set calcu- 
lation. The total energy used in DFT calculations is 

E[n] = Ts[n] + j U„„(r)n(r)dr + i;„„({R,}) 

+ iy" n{v)n{r')V{Yy)dvdv' + E,,[n]. (1) 

The kinetic energy Tg [n] and exchange-correlation en- 
ergy Exc[n] are in essence local (i.e., do not interact 
across cells) and therefore do not need to be changed.El 
The terms through which the contents of different super- 
cells interact are the Hartree energy, the ionic potential 
Vion and the ion-ion Coulomb repulsion energy Eion- In 
an infinite crystal these three terms give rise at-jjecipro- 
cal lattice vector G = to infinite contributionfjl^ which 
can be shown to cancel. 

We first consider the electron-ion term. The ionic po- 
tential felt by the electrons is given by 



l^.o„(r)=^^U„(r + l) 



(2) 



where the sums are over all lattice vectors 1 in the infinite 
crystal and ions a in the unit cell, and Va is the pseudopo- 
tential for atom a. In order to carry out a calculation 
with the Coulomb interaction cut off on the supercell, we 
simply require the sum over the contents of the unit cell 
only. Whilst in the infinite case the sum over unit cells 
means that we require the local pseudopotential only at 



the reciprocal lattice vectors (all other components in the 
sum cancel), we now need to integrate the local poten- 
tial in reciprocal space up to the energy cutoff and then 
FFT to the correct real-space ionic potential for the finite 
cell. We use a localized real-space representation of the 
non-local part of the pseudopotential,lliJ which is readily 
modified to prevent any interaction across cells. 

The ion-ion energy is now evaluated easily by a direct 
sum over the ions_in one cell (as opposed to the usual 
Ewald summationEj). 

The remaining term through which supercells interact 
is the Hartree potential, given by 



VH{r) 



n(r') 



dr'. 



(3) 



all space 



As this is an infinite convolution it is normally evaluated 
in reciprocal space where 



Vh{G) = n(G)U(G) 



(4) 



where V^(G) is the Fourier transform of the Coulomb 
interaction, U(G) = 47r/|G|2. 

To remove electrostatic interaction between electronic 
charge densities in different supercells we require the 
modified Hartree potential, 



(r) = / ■ 

JSC 



(5) 



where the integral is performed only over the cell con- 
taining r. This effectively defines a modified Coulomb 
interaction. 




for r, r' in same cell 
otherwise. 



(6) 



This cutoff can be imposed either in real or reciprocal 
space, and we test the effect of each in the next sec- 
tion. Such an interaction is clearly no longer a function 
of r — r' only though, and hence does not preserve the 
simple multiplicative form of (^). Onida et al. employed 
a spherically cutoff Coulomb interaction. 
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U(r,r') = 



|r - r' 




for |r - r'l < i? 

otherwise 



(7) 



for GW calculationfE on negatively charged clusters to 
prevent screening by periodic images in other supercells. 
As noted in their work, though, an interaction cut off 
only as a function of |r — r'| clearly in principle requires 
additional separation of densities in adjacent supercells 
to remove all interaction effects. 

An interaction with the form of ^ can still be used 
in conjunction with the convolution theorem, by pacLding 
the density with an external region of zero density.L3 If 
we consider, for simplicity, a cubic supercell of side L, 
then we take a new grid of side 2L, defining the density 
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to be zero on the extra points and evaluating the usual 
real-space Coulomb interaction on the expanded grid. It 
can then be seen that performing an FFT on the padded 
density and interaction, multiplying them as in and 
reverse-transforming to the padded cell allows us to ex- 
tract the correct Hartree potential [as defined in (||)] in 
the original real-space cell. For the numerical transform 
of the Coulomb interaction the singularitii-|at r = r' is 
integrated and averaged over one grid unitJlil A spherical 
cutoff in the Coulomb interaction at a range of -^/S L (to 
span the cell) facilitates an analytical transform to recip- 
rocal space but also means that the padded grid must 
now have side (1 -I- -^3) L to strictly avoid any spurious 
contributions to the Hartree integral. 

This approach to cutting off supercell interaction in a 
plane-wave calculation thus leads to a procedure for per- 
forming the Hartree integral identical to that employed 
in some finite-grid calculations, e.g., for molecular dy- 
namics of charged clustersJlj Viewing this procedure as 
a physical truncation of the Coulomb interaction in an 
infinite calculation, rather than as a mathematical tool 
to perform a finite integral, makes the potential for sav- 
ings in computational effort more apparent. In Fig. |l| we 
show electron density in a unit cell of side L padded by 
a distance D and interacting with all density within a 
box of side 2R. If D and R are set equal to L then we 
obtain the correct answer for the finite Hartree energy. If 
R is reduced this corresponds to cutting off interaction 
between electronic density on opposite sides of the same 
cell. D can always be reduced to R, but if it is reduced 
still further then spurious interaction between density on 
the edges of adjacent supercells is introduced. As the 
density will be exponentially small at the box edges for a 
converged calculation, it is conceivable that reducing R, 
or D below R, will only introduce very small negative 
and positive errors respectively in the total energy. As 
this would be a useful procedure in cases where calcula- 
tion of the Hartree energy dominates computation time, 
we test the extent to which such savings can be made in 
the next section. 

Through these relatively straightforward alterations to 
the standard infinite code, a finite calculation can be per- 
formed with all electrostatic interaction between super- 
cells removed, and with all other aspects of the method 
(pseudopotentials, conjugate gradients, geometry optimi- 
sation) unchanged. 

III. RESULTS 

A. Test systems 

We first examined the test case of a charged system, 
taking the ionisation potential of a Mg atom as used 
by Makov and Payne to assess the convergence of their 
energy-correction method. The cutoff Coulomb interac- 
tion allows a charged system to be treated in exactly 



the same way as a neutral system. For a usual supercell 
calculation it is necessary to insert a neutralising back- 
ground to remove the divergent repulsion energy of the 
charged lattice of cells. The binding energy between this 
background and a lattice of point charges is then ac- 
counted for by the Madelung correction, whilst Makov 
and Payne gave an expression for a further correction in 
terms of the charge density profile. 

For the neutral atom, use of the cutoff interaction gave 
almost identical (difference of order 10~^ eV) results to 
the unmoderated supercell calculation for all box sizes 
above 8 A as expected: the Hartree, electron-ion and 
ion-ion energy terms for the finite and repeated supercell 
systems sum to the same total in this case, as the atom 
has no non-zero moments. 

The calculation for the positive ion with the cutoff in- 
teraction can be undertaken simply by reducing the elec- 
tron number by one as all interaction between supercells 
has been removed. This gave an energy for the ion which 
was converged within 0.01 eV at a box size of 8 A, as op- 
posed to 10 A for equivalent convergence of the energy- 
correction method. We note that there is no power-law 
behaviour in the convergence of our calculation, and the 
answer is in principle exactly that of the true finite sys- 
tem once the box is large enough for the electron density 
to be zero at the edges. 

We next consider a system with a dipolc moment. In 
Fig. H we show the potential energy of stretching a NaCl 
molecule by 0.3 A from its calculated equilibrium length 
of 2.231 A, as a function of cubic box size used for the 
calculation. The removal of all electrostatic interaction 
between supercells eliminates the slow power-law conver- 
gence, and also represents some improvement over the 
convergence with respect to supercell size achieved by 
the energy-correction approach. 

We also tested the effect of reducing the interaction 
range R and the padding distance D as outlined in the 
previous section. At a cubic box size of side 12 A the 
total energy for the Mg atom is converged to within a 
few meV. We found that the range parameter could be 
reduced below 0.6L, and the padding distance then fur- 
ther reduced below 0.31/, without either change affecting 
the total energy by more than 10"'' eV. As this gives 
a padded grid of side 1.31/ (as opposed to 2L in the 
exact unaliased case) this can actually give an improve- 
ment of almost an order of magnitude in the time taken 
to compute the Hartree energy. We therefore make use 
of reduced range and padding for all further calculations, 
whilst restricting the resultant error in the total energies 
to be less than 10"'^ eV. 



B. Polar molecules 

The test cases of a charged and dipolar system con- 
firm that the cutoff interaction removes power-law con- 
vergence and indicate that it can yield some savings com- 
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pared to finite-grid calculations. However, in the case of 
the NaCl bond, the energy difference being measured is 
rather small (around 0.2 eV), while the system itself con- 
sists of an extremely large dipole (calculated value of 8.2 
D) in a fairly small box. As an example of a more realis- 
tic system of interest we calculated the structure of the 
molecule dimethylsulphoxide (DMSO), shown schemati- 
cally in Fig. H, which still has a relatively large dipole 
(around 3.9 D), using the normal supercell approach and 
the cutoff interaction. Both calculations were converged 
to within 0.01 eV per atom, and to within 0.01 eV per 
atom of each other, for a cubic box size of 8 A and an 
energy cutoff of 650 eV. As shown in Table ||, by the time 
the calculation is converged to this degree, the structures 
obtained already agree to within 0.002 A for bond lengths 
and 0.3° for bond angles. This is a surprising result given 
the large dipole moment of DMSO, and the relatively 
small distance by which the molecules in the supercell 
calculation are separated. 

There is also some interest in thajcalculation of molec- 
ular dipole moments using DFT.EEl We found that the 
calculated dipole showed a greater proportional differ- 
ence between methods than the structure, as a result of 
neither calculation giving a particularly well-converged 
dipole for box sizes of less than 12 A, where the dipoles 
agree within 0.1 D. This is because the dipole is far more 
sensitive to small movements of charge density from one 
side of the supercell to the other than the total energy. 

Use of the fully unaliased cutoff interaction made the 
calculation around five times slower than the unmoder- 
ated supercell calculation for DMSO in the 8 A box. The 
extra computation time is used almost entirely in the 
EFT of the padded density. Fully exploiting the reduced 
padding and range, by treating them as additional con- 
vergence parameters, as discussed earlier, resulted in a 
cutoff calculation now approximately twice as slow. 

The cutoff interaction method (or a finite-grid calcu- 
lation) does not rely on any particular choice of super- 
cell geometry. In the case of an infinite cubic lattice of 
point dipoles, the dipole-dipole interaction automatically 
cancels. Makov and Payne derived an energy correction 
for the macroscopic field imposed by PBCs and thus ob- 
tained the energy of the infinite cubic array which, due 
to the dipole-dipole cancellation, converged to the re- 
sult for the isolated system as order . For non-cubic 
boxes the dipole-dipole interaction does not cancel and 
larger errors in a supercell calculation might be expected. 
The use of an artificially cubic supercell for elongated 
molecules in a supercell calculation, if needed, would in- 
volve expensive additional vacuum. 

We therefore next considered the case of an elongated 
molecule in a non-cubic supercell. We performed a ge- 
ometry optimisation for the molecule p-nitroanisole, with 
structure as shown in Fig. ^ The box size for our con- 
verged calculation was 14 A x 10 A x 6 A, and an energy 
cutoff of 650 eV was needed. Representative results ob- 
tained for bond lengths and angles are shown in Table ||. 
Once again the similarity is striking, with bond lengths 



agreeing within 0.001 A and bond angles within 0.1°. 
Thus the structure of p-nitroanisole is described equally 
well by the supercell technique despite the non-cubic ge- 
ometry (and even though the molecule has a somewhat 
larger dipole, calculated value 5.2 D, than DMSO). Total 
energies with and without interaction between molecules 
also agreed to well within 0.01 eV per atom. We empha- 
size that (even for these organic molecules chosen specif- 
ically for their atypically large dipole moments), in or- 
der to obtain this agreement between supercell and finite 
results, the supercell was only as large as was required 
for convergence of the energy of the isolated system i.e., 
large enough to accommodate the charge density of the 
molecule. 



C. Charged clusters 

As outlined earlier, supercell calculations for charged 
systems require the energy to be corrected by the in- 
sertion of a neutralising background, the Madelung cor- 
rection, and, if necessary, the Makov-Payne correction, 
depending on the degree of accuracy required. The first 
two terms affect only the G = component of the to- 
tal potential. This means that an unmoderated supercell 
calculation can in fact be carried out for a structure de- 
termination (because for a fixed number of electrons the 
error in total energy due to the first two terms is geome- 
try independent). Any remaining errors in the resulting 
structure would be caused by interaction effects of the 
order addressed by the Makov-Payne correction. 

We therefore carried out a structure determination for 
the Na2^ cluster. We obtained a converged bond length 
with the cutoff interaction of 3.513 A, and with the su- 
percell method of 3.512 A, in a box with dimensions 
12Axl0Axl0A, the smallest box size at which either 
calculation was well converged. In this case also, interac- 
tion effects (beyond the straightforward interaction of a 
lattice of charges which can be easily removed) have very 
little effect on the structure determination. 



D. Coulomb interaction in finite-grid calculations 

In a PBCs calculation the Coulomb interaction is nec- 
essarily defined in reciprocal space - although the interac- 
tion is infinitely ranged both in real and reciprocal space, 
components of the density only extend out to some cutoff 
in reciprocal space, allowing in principle an exact deter- 
mination of the Hartree energy. 

In a finite-grid calculation the Coulomb interaction, 
which is now only required over some limited real space 
volume, is usually defined in real space as 1/r. However, 
when viewed as a cutoff interaction in a periodic system, 
it is possible to use forms specified analytically in either 
real or reciprocal space. 
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As an example, we performed calculations for a Mg 
atom, using a Coulomb interaction cut off spherically at 
radius R = y/3L in real space, with appropriate padding, 
to determine the Hartree energy. This cutoff was im- 
posed first by taking the analytical Fourier transform to 
reciprocal space. 



V{G) 



1 - cos{\G\R) 



(8) 



and then in real space, taking 1/r within the sphere and 
zero outside. Clearly the FFT of each should give the 
other, but only in the limit of a large energy cutoff (i.e., 
fine Fourier transform grid). 

Convergence of the resulting total energies is shown in 
Fig. p. ft can be seen that in fact convergence is much im- 
proved by using a cutoff interaction specified analytically 
in reciprocal space, rather than using 1/r with averaging 
of the singularity at the origin. 



IV. DISCUSSION 

Our calculations using the cutoff Coulomb interaction 
allow for detailed and unambiguous assessment of the ef- 
fects caused by interaction between supercells. As the 
errors in energy caused by such interaction effects are 
relatively small, it is a crucial aspect of our work that we 
are able to retain the same pseudopotentials, exchange- 
correlation approximation, and basis set for the super- 
cell and cutoff calculations. Although image interaction 
has often been cited as an argument against the super- 
cell method, these factors have generally prevented accu- 
rate or conclusive comparisons, and the ability to perform 
such a study is the main motivation for our approach. 

We firstly make two points regarding the implications 
of our results with the cutoff Hamiltonian for finite-grid 
calculations. The savings we obtain through the reduced 
interaction range and padding verify that, although it is 
important to converge the exponential tail of the den- 
sity for the total energy, the contribution to the Hartree 
energy from the interaction of these tails across one cell 
or between adjacent cells is not as significant _ This also 
explains why the calculations of Onida et aln obtained 
well-converged results without employing padding, and 
without needing to fully separate clusters in different unit 
cells. This procedure is certainly worth considering in 
finite-grid calculations, although the proportional saving 
will be less for large systems (where the density decays 
exponentially only over a smaller proportion of the cell) . 

Implementation of the Coulomb cutoff in real and re- 
ciprocal space showed a marked difference in convergence 
with respect to energy cutoff. This is because the recip- 
rocal space form treats the singularity in 1/r analytically. 
The slow convergence seen in the total energy with the 
real space interaction is simply a result of the need to 
sample 1/r very finely, and does not relate to the energy 
cutoff required to converge the density, which is far lower. 



As a result a cutoff interaction defined in reciprocal space 
(and then transformed to real space if necessary) should 
offer improved convergence for any finite computation of 
the Hartree energy. 

The results of our calculations, especially for the dipo- 
lar molecules, reveal a very high level of accuracy in 
the structures and energies yielded by the supercell 
method. In particular, the common assertion that ad- 
ditional vacuum is required proves not to be the case for 
the charged cluster Na2^ or for the molecules DMSO 
and p-nitroanisole, even though both have large dipoles. 
The excellent agreement seen in energies and structures 
suggests that the additional effort involved in a fully fi- 
nite calculation is unlikely to be necessary for a wide 
range of clusters and molecules, as these errors are con- 
siderably smaller than the differences given by choice of 
pseudopotential or exchange-correlation approximation. 
This conclusion is consistent with the supercell calcula- 
tions on Cu~ clusters by Massobrio et alE!\ who, using a 
positive jellium background, found that the size of the re- 
maining error caused by interaction, as evaluated by the 
Makov-Payne formula, was already negligible compared 
to total-energy differences. 

We therefore suggest that the most efficient method 
for determining structures and energies of clusters and 
molecules is to perform a conventional band structure 
calculation, converged with respect to the supercell di- 
mensions. For charged systems or dipoles on a cubic lat- 
tice, involving, for example, very small energy differences, 
or very large dipoles, corrections to the energy can be as- 
sessed using the Makov-Payne terms. We note, though, 
that this approach for dipolar molecules can only be used 
in conjunction with cubic supercells. Another option is 
to restart from the optimized geometry with the cutoff 
interaction. We re-iterate though that, as we have found 
errors related to interaction to be less than 0.01 eV per 
atom in realistic cases, we expect that an uncorrected 
supercell calculation should nearly always suffice. 

Most finite calculations are undertaken with Gaus- 
sian basis sets. In the absence of any adverse effects 
from enforced periodicity, plane-wave basis sets offer 
some advantages for ah initio calculations of clusters and 
molecules. Results of calculations undertaken with Gaus- 
sian basis sets are often dependent on the choice of basis 
set, as a result of being difficult to converge with re- 
spect to basis set size: there are also problems associ- 
ated with lack of basis set orthonormality and descrip- 
tion of diffuse occupied states. Plane waves offer a simple 
and systematic convergence procedure, and also facili- 
tate straightforward computation of forces for molecu- 
lar dynamics, and there has therefore been considerable 
interest in the applicability of plane-wave basis spts for 
first-xirinciples chemical calculations for moleculeala and 
ions.Ej Our work provides clear evidence that interaction 
between supercells provides no obstacles to a plane- wave 
approach. 

Interaction effects are also potentially an issue in calcu- 
lations of infinite systems such as polymers or surfaces. 
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where a finite-grid calculation is not possible. It often 
proves necessary to include several la yecs .of vacuum in 
supercell slab calculations of surfaces jESEil for example. 
In the light of our results it seems possible that this may 
be purely the amount of vacuum needed to describe the 
exponential tail of the density, rather than an electro- 
static effect. This would mean that no larger a supercell 
would be needed for the supercell method than would be 
required by using localized orbitals, although the differ- 
ing geometry and larger polarizability do allow for po- 
tentially larger (and physically different) errors in these 
cases. 

Similarly, in order to simulate defects in solids a very 
large supercell may be needed. Makov and Payne de- 
veloped an energy correction fot.4his case, based on ear- 
lier work by Leslie and GillanJ23 but our findings also 
raise the question of whether differences from the answer 
for the single defect in an infinite crystal are dominated 
by the need to describe the physical environment of the 
defect over medium range length scales or indirect elec- 
trostatic effects, rather than by direct electrostatic in- 
teraction between defects. Poor convergence because of 
chemical interactions is consistent with the finding that 
convergence wi1|t, respect to supercell size depends on k- 
point sampling.E2l 

The Makov-Payne correction to the total energy for 
dipolar systems in the case of cubic supercells accounts 
for much of the (small) difference between the PBCs cal- 
culation and the cutoff interaction, as was shown in Fig. 
||. The correction has the effect of reinstating the macro- 
scopic field of the infinite array of supercells (which has 
been removed by the depolarizing field implicit in the 
use of PBCs), and calculates the corresponding energy. 
The remaining difference between the corrected energy 
and the isolated system is the interaction energy of the 
lattice. As the energy of a cubic lattice of point dipoles 
is zero, the corrected energy per cell should then give a 
good approximation to the energy of the finite system. 
We verified explicitly that the remaining (very small) 
discrepancy between the Makov-Payne-corrected energy 
and the converged result for the finite system could be 
accounted for almost entirely by the energy of a classical 
lattice of appropriate finite-length dipoles. It is also use- 
ful to observe that, although strictly the Makov-Payne 
correction requires a modification of the electrostatic en- 
ergy functional and accordingly the effective potential, 
simply applying the correction using the density from 
a conventional supercell calculation actually provides a 
good estimate of the total-energy error. 

For an infinite non-cubic lattice of dipoles, we cannot 
generally expect the interaction energy to be zero. How- 
ever, for the p-nitroanisole molecule we studied, the in- 
teraction energy of the lattice (again evaluated for finite 
length dipoles) can be seen to have the opposite sign of 
the correction for the depolarising field, thus explaining 
why the supercell calculation for this system also differed 
only slightly from the cutoff calculation. 



V. CONCLUSIONS 

We have used a modified Coulomb interaction, cut off 
in real space, to allow calculations for finite systems to 
be undertaken within periodic boundary conditions, with 
no interaction between supercells. We have demonstrated 
that, in any cases where a calculation for an isolated sys- 
tem should prove necessary, the cutoff Coulomb interac- 
tion technique can provide considerable efficiency gains 
compared with a finite-grid calculation by allowing a de- 
gree of aliasing in the Hartree integral without affecting 
the overall answer. We have also found that convergence 
with respect to energy cutoff of finite-grid calculations 
can be improved by defining the Coulomb interaction an- 
alytically in reciprocal space. 

Detailed comparison of calculations performed for 
charged clusters and molecules with large dipoles shows 
that the supercell method gives very accurate energies 
and structures compared with the isolated system, mean- 
ing that the extra expense of finite-grid calculations is 
unnecessary in the large majority of cases. This shows 
that plane- wave basis set calculations are applicable to a 
wide range of chemical and biochemical problems with- 
out the need for artificially large supercells or corrections 
for supercell interaction. 
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L (Angstroms) 

FIG. 2. Energy of a NaCl bond with an extended bond 
length relative to equilibrium as a function of box size L, 
from suporccll method (circles), corrected by Makov-Paync 
functional (diamonds), and using the cutoff Coulomb interac- 
tion (triangles). 
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FIG. 3. Structure of DMSO, which has a dipole as shown. 
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FIG. 1. Schematic illustration of padding procedure for cal- 
culating the Hartree energy. Electronic density in a unit cell 
of side L is padded by a distance D and interacts over a box 
with range parameter R (see text). 



FIG. 4. Structure of p-nitroanisole, with dipole as indicated. 
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FIG. 5. Total energy of Mg atom, calculated with cutoff 
Coulomb interaction, defined in real space (triangles) and re- 
ciprocal space (circles). 
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TABLE I. Bond lengths and bond angles for the molecule 
DMSO in a cubic box of side 8 A calculated using the usual 
supercell technique and with cutoff electrostatic interaction. 





S-O 


S-C 


C-S-C 


C-S-0 


Supercell 


1.433 A 


1.785 A 


95.8° 


106.7° 


Cutoff 


1.431 A 


1.786 A 


96.1° 


106.6° 



TABLE II. Bond lengths and bond angles for the molecule 
p-nitroanisole in a 14 A x 10 A x 6 A box calculated using the 
usual supercell technique and with cutoff electrostatic inter- 
action. 





N-C 


O-C 


N-0 


O-N-0 


Supercell 


1.485 A 


1.267 A 


1.078 A 


125.7° 


Cutoff 


1.485 A 


1.266 A 


1.078 A 


125.7° 
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